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ABSTRACT 

I introduce batman, a Python package for modeling exoplanet transit and 
eclipse light curves. The batman package supports calculation of light curves 
for any radially symmetric stellar limb darkening law, using a new integration 
algorithm for models that cannot be quickly calculated analytically. The code 
uses C extension modules to speed up model calculation and is parallelized with 
OpenMP. For a typical light curve with 100 data points in transit, batman can 
calculate one million quadratic limb-darkened models in 30 seconds with a sin¬ 
gle 1.7 GHz Intel Core i5 processor. The same calculation takes seven minutes 
using the four-parameter nonlinear limb darkening model (computed to 1 ppm 
accuracy). Maximum truncation error for integrated models is an input pa¬ 
rameter that can be set as low as 0.001 ppm, ensuring that the community is 
prepared for the precise transit light curves we anticipate measuring with up¬ 
coming facilities. The batman package is open source and publicly available at 
https://github.com/lkreidberg/batman, 

Subject headings: methods: data analysis - methods: numerical 


Introduction 


The transit technique has revolutionized the study of exoplanetary systems. Thanks 
l argely to the Kep ler mission, thousands of planets have been discovered with this method 
flRowe et al.l 120151) . These discoveries have yielded transformative constraints on planet 
occu rrence rates over a wide range of planet sizes, orbital periods, and host star proper- 
ties f Borucki et ah 2011 : Youdin 2011 : Howard et ah 2012 : Fressin et ah 2013 : Dong fc Zhu 
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2013 
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2013; 

KoDDaranu 

2013; 

Foreman-Mackev et al. 

2014; 

Dressing & Gharbonneau 

201 bh. Tran- 


sit light cu rves can even reveal pl a nets’ atmospheric temperature structure and compo- 


sition 

(e.g. 

Seager & Sasselov 

2 OO 0 I: 

Gharbonneau et al.l 

2002: 

Lecavelier Des Etangs et al. 

2008 


Sing et a’ 

2011 


Deming et al. 

2013 


Knutson et al. 

2014a 


Fraine et al. 

2014; 

Kreidberg et al. 

2015 

). A number of current and planned observational facilities - including 


K2, TESS, CHEOPS, JWST, and PLATO - will measure precise transit light curves for 
thousands of exoplanets that will further advance our understanding of planet formation, 
evolution, and habitability. 

Light curve models are a fundamental tool for transiting exoplanet science, but they are 
not trivial to compute quickly and accurately. Accurate calculation is challenging because 
the model must account for the planet’s size and position on the sky, as well as stellar limb 
darkening, which causes the apparent brightness of the stellar disk to decrease from center to 
edge. The stellar inte nsity prohle can be ht with sev eral function al forms, incl uding a linear 

sqnare-root 
), exponential 


l imb darkening law fISchwarzschild fc Villigerl Il906l). quadratic flKopal 


f Piaz-Cordoves fc Gimenezlll992l) . logarithmic flKlinglesmith fc Sobieski 


( Claret fc Hauschildt 20031) . and four-parameter nonlinear fjciaretl 200oh. For some of these 


1951 


197C 


prohles, mode l transit light curves can be calculated analytically flMandel fc Agoll 12002 


GimenezI I 2 OO 6 I: lAbubekerov fc Gostev! l2013l) . Other prohles do not have analytic solutions, 
and models must be calculated by numeric integration of the stellar intensity over the disk 
of the planet. In addition, speed is an important consideration because a large number of 
models must typically be calculated to make a robust estimation of transit parameters and 
their uncertainties. 


A number of codes are available to calculate transit light curves. iMandel fc Agoll (120021) 
provide Fortran and IDL rontines to comp ute models for qua dratic and nonli near limb dark¬ 
ening laws. The software packages TAP (IGrazak et al.N201^ and EXOFAST (lEastman et al. 


20131) include IDL implementations of the iMandel fc Agoll (120021) algorithm for quadratic 
limb darkening. JKTEB(1P calc ul ates models in For t ran fo r a broad range of limb darkening 
laws (jSouthworth et al.l 120041) . iKiurkchieva et aLl (l2013l) introduce the pure Python code 
TAC-maker, which performs numeric integration for arbitrary limb darkening prohle s. There 
are a l so routine s available to m odel simult a neous transits by one or more bodies (IKipping 
2 OIII: |Palll20121) . Most recently. iParviainenI (l2015h releas e d the Pyt hon package PyT ransit, 
which implements analytic models from IMandel fc Agoll (120021) and I GimenezI (120061) . 


In this paper, I introduce the open source Python package batman. This package is 
based on code that was used to model high-precision light curves obtained for atmosphere 
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characterization ( Kreidberg et al. 2014alfl boiS : Stevenson et al. 2014a b 313). The batman 
package enables fast computation of transit light curves for any radially symmetric limb 
darkening law, and currently supports uniform, linear, quadratic, logarithmic, exponential, 
and nonlinear limb darkening. Ligh t curves for the hrst th ree of these are calculated analyt¬ 
ically based on the formalism from iMandel fc Agoll fj2002[) . Models for the remaining cases 
are computed with an efficient new integration scheme, described in §|2l The package also 
supports secondary eclipse modeling. I discuss batman’s features and performance in §[3] and 
conclude in §01 


2. Algorithm 

To calculate the fraction 6 of stellar flux blocked by a transiting planet, one must 
integrate the sky-projected intensity of the star (J) over the area obscured by the disk of the 
planet (S'): 

5 = IJ MS (1) 

where I is normalized such that the integrated intensity over the stellar disk is unity. This 
expression is valid for any general stellar surface brightness map; however, it is slow to 
evaluate numerically because the differential area elements must be small (< 10“®) in order 
to achieve better than one part per million (ppm) accuracy. 

On the other hand, if the stellar intensity prohle is radially symmetric, the two- 
dimensional calculation in Equation[T] can be reduced to one dimension and sped up greatly 
with the following algorithm: 

5 = '^ I [A{xi,rp,d) - A{xi_i,rp, d)] (2) 

i=l ^ ' 

where x is the normalized radial coordinate 0 < a; < 1, I{x) is the ID stellar intensity prohle, 
rp is the planetary radius (in units of stellar radii), d is the separation of centers between the 
star and the planet (in stellar radii), and A{x,rp,d) is the area of intersection between two 
circles of radii x and r^, separated by a distance d. The sum is carried out over the range 
Xq = MAX((i — Tp, 0) to Xn = MIN((i -|- r^, 1). The intersecting area is given by: 


COS“^ U + Tp COS“^ V — 0.5y/w, 

Tp — d<x<rp + d 


TTX^, 

X < Tp — d 

(3) 

Tir^ 

P 

X >rp + d 
























where 
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( 4 ) 

( 5 ) 

( 6 ) 


u = {d^ + — Tp) / {2dx) 

V = + Tp - x^)/{2drp) 

w = {—d + x + rp){d + x — rp){d — x + rp){d + x + rp). 


See Figured] for a schematic illustrating the geometry of the integration. The advantage 
of this integration scheme is that the only error introduced is due to approximating the stellar 
intensity as constant over the differential area element AA = A{xi,rp,d) — A{xi-i,rp,d). 
Computation of AA is expensive, but it needs to be calculated relatively few times (<< 10® 
for sub-ppm accuracy). This makes the integration faster than a scheme with a simpler area 
element (e.g., A A = Ax Ay), that requires a much smaller step size to achieve the same 
accuracy. 


The integration can be further optimized by using a nonuniform step size. Typical stellar 


inten sity prohles have larger gradients near the limb of the star than at the center (e.g. I Claret 


20001 ). so smaller steps are required at larger x values to achieve the same accuracy. I adopt 


the following step-size scaling: 


Xi - Xi-i = / cos ^ {Xi-i) 

where / is a constant scale factor. This prescription is fast to compute and well-behaved at 
the limits a; = 0 and a; = 1. 


3. The batman package 

The Python package batman implements the algorithm described in §[2] and sev¬ 
eral analytic models to calculate transit light curves. batman is an open source 
project and is being developed on GitHub. Full documentation is available at 
https://github.com/lkreidberg/batman. I summarize the main capabilities of the pack¬ 
age here. 


3.1. Limb Darkening Models 


batman supports calculation of exoplanet transit light curves for uniform, linear, 
quadratic, square-root, logarithmic, exponential, and four-parameter nonlinear stellar in- 
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Fig. 1.— Schematic illustration of the integration scheme. The star (large black circle; 
partially visible) has a radius of 1 and is centered in the plane of the sky at (x, y) = (0, 
0). The planet (smaller black circle) is separated from the center of the stellar disk by a 
distance d (marked by the solid black line). The star is partitioned into concentric circles 
(dotted lines) in order to calculate the integral over the planet disk. A single integration 
element AA is shaded in orange. The integration step size illustrated here is larger than 
for a typical calculation for visual clarity. Note that because the stellar intensity prohle is 
radially symmetric, the coordinate system can be chosen snch that the planet lies on the 
x-axis, as shown. 
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tensity profiles: 


J(/i) = Jo 

(uniform) 

(7) 

/(/i) = /o[l-Ci(l-/i)] 

(linear) 

( 8 ) 

/(/i) = Jo[l - Ci(l - yu) - C2(l - nY] 

(quadratic) 

(9) 

/(/i) = /o[l - Ci(l - yU) - C2(l - y/jl)] 

(square-root) 

(10) 

/(yfi) = /o[l - Ci(l - yU) - C2yUlnyU] 

(logarithmic) 

( 11 ) 

/(yfi) = Jo [1 - ci(l - yu) - C 2 /(l - expyu)] 

(exponential) 

(12) 

/(yfi) = Jo[l - Ci(l - yfi^/^) - C2(l - yfi) - C3(l - yfi^/^) - C4(l - 

yfi^)] (nonlinear) 

(13) 

where /i = y/l — and ci,..., c„ are limb darkening coefficients. 

The batman source 

distri- 


bution also includes a template for the creation of a custom profile for any radially symmetric 
function. 


The square-root, logarithmic, exponential, nonlinear, and custom models are computed 
with the numeric integration scheme from §[2l The uniform, linear, and quadratic models 
are calculated analytically, wit h code based o n the Fortran routines occultquad.f and 
occultuniform. f pr ovided by Mandel fc Agol f 200211 . For the a nalytic models, I follow 


Eastman et ahl (120131) and use the algorithm from iRnlirschl (119651) to improve calculation 


speed and accuracy for elliptic integrals of the third kind. 


3.2. Secondary Eclipse Model 

batman can also model secondary eclipses. Eclipse light curves are generated with 

f = 1 + fp{l- a) 

where / is normalized flux, fp is the planet-to-star flux ratio, and a is the fraction of the 
planet disk that is occulted by the star. The model is normalized such that the stellar flux 
is unity. For a separation d, the occultation fraction a{d) = at{d)/rp, where 1 — at{d) is the 
transit light curve with uniform limb darkening. Note that this model assumes the planet 
flux is constant for all orbital phases. 


3.3. Utilities 

batman includes a utility function to calculate the separation of centers d between the 
star and the planet based on orbital parameters of the system. The input parameters are 
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the planet semi-major axis a, inclination i, eccentricity e, longitude of periastron ca, period 
P, and time of inferior conjunction to. The separation of centers is given by: 


d = 


a{l 


1 — sin^ (cj -I- /) sin^ i 


1 -|- e cos / 

where / is the true anomaly. The true anomaly is calculated with the algorithm provided by 
Murray & Correia in Chapter 1 of lSeagerl (120101) . For circular orbits, we adopt the convention 
/ = 7r/2. 


batman provides utilities to calculate the time of periastron, time of inferior conjunction, 
and the se condary eclipse t ime f rom the other transit parameters, using the method described 
in §3.1 of lEastman et al.l (120131) . Note however that batman does not correct for the effects 
of light travel time. 


An additional utility is light curve supersampling, which allows the user to calculate the 
average value of the light curve over a specihed number of evenly spaced points during an 
exposure. 


3.4. Accuracy 


Recent transit observat i ons have yielded signa l-to-noise greater than 1000 per exposure 


(e.g. iKreidberg et al.ll2014bl: Knutson et al.ll2014bl) . Accurate transit light curve calculation 
is essential for modeling such high precision measurements and will be increasingly important 
for data obtained with next-generation facilities, batman therefore enables the user to specify 
the maximum allowable truncation error for numeric integration. Figure[2] shows an example 
transit light curve and its truncation error. 


To ensure that the truncation error is below the specihed threshold, the integration 
step size is tuned during model initialization. Truncation error is measured relative to a 
model calculated with a very small step size (/ = 5 x 10“'^). For typical limb darkening 
prohles, this method is reliable for truncation errors down to ~ 10“^ ppm. However, tuning 
the step size is a slow operation because it requires computing several light curve models 
(~ 10). As an alternative, methods are available to set the step size directly and calculate 
the corresponding truncation error. 

I tested the accuracy of the analytic model for quadratic limb darkening by comparing 
it to a numerically integrated model with an error tolerance of 0.001 ppm. The analytic 
model is accurate to 0.03 ppm for a test case with = 0.1, (ci,C 2 ) = (0.1, 0.3), sampled at 
10® evenly spaced points over the interval 0 < d < 1. The accuracy is somewhat worse than 
machine epsilon because of error tolerance in the computation of special functions. 

















I also tested the accuracy of the widely-used iMandel fc Agoll (1200211 code that uses 
Numerical Recipes functions to calculate elliptic integrals flPress et al.lll992[) . The accuracy 
was better than 0.005 ppm for most input values; how ever, for t he ca se Vp — d < e, the error 
in the light curve exceeded 2 ppm. By contrast, the iBulirschl fjl965[) algorithm for elliptic 
integrals is well-behaved for this case and also faster. 


3.5. Performance 


Computationally intensive sections of code (including all of the transit model calcu¬ 
lation) are written as C extension modules with the Python/C API, which improves the 
performance by a factor of 30 over a pure Python implementation for quadratic limb dark¬ 
ening. batman also includes the option to parallelize at the C level with OpenMP, which 
further speeds up the calculation. The number of processors is specihed by the user, batman 
will raise an exception if the user attempts to parallelize a calculation on a system where 
OpenMP is not supported. 


I tested batman’s performance over a range of typical use cases with a 1.7 GHz Intel 
Core i5 processor. In Figure[3l I show the truncation error versus execution time for a 
single transit light curve calculation using a nonlinear intensity prohle, compared to the 
execution time for a quadratic model computed analytically. The test case consists of 100 
points evenly sampled in time during the planet’s tran sit. I used physical parameters for 
the transiting planet GJ 1214b (iKreidberg et al.ll2014bl) . The stellar intensity prohle is the 
same for the nonlinear and quadratic models: the nonlinear limb darkening coefficients are 
(0.0, 0.7, 0.0,—0.3) and the quadratic coefficients are (0.1, 0.3). Decreasing the truncation 
error by a factor of 10 increases the computation time by a factor of three. 


3.6. Comparison with Analytic Models for Nonlinear Limb Darkening 


I explored using an analytic model to calculate transit light curves for the four-parameter 
nonlinear lin ib darkening pr o hle. The analytic solution for nonlinear limb darkening was 
presented in iMandel fc Agoll ( 20021). but it is not us ed in any published software packages. 
The original code provided bv IMandel fc Agoll (120021) uses a numeric integration scheme that 
is 20 times slower than the algorithm presented in §[2] (for an error tolerance of 1 ppm). 


The analytic solution is challenging to compute because it uses the Appell FI hyper¬ 
geometric function. This function is only convergent for certain reg ions of parameter space 
and must be calculated with analytic continuation for other cases. IColavecchia fc Gasaneo 
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Time from central transit (hours) 


Fig. 2.— An example transit light curve for a nonlinear stellar intensity profile (top panel) 
and truncation error for the calculation (bottom panel). The error tolerance parameter was 
set to 1.0 ppm. The truncation error increases with distance from the center of the star up 
to around ±0.2 hours from the time of mid-transit, because the stellar intensity gradient 
is larger at larger radii. The error decreases again during ingress and egress as the planet 
eclipses a smaller fraction of the stellar disk. 
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Execution time (s) 


Fig. 3.— Truncation error as a function of execution time for a light curve modeled with the 
nonlinear limb darkening law (black line). The execution time for a quadratic model (com¬ 
puted analytically to better than 0.05 ppm accuracy) is indicated by the arrow. Calculations 
were made with a 1.7 GHz Intel Core i5 processor. 
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fl2004j) provide a Fortran library for computing Appell FI. I used this library to implement 
the analytic model for nonlinear limb darkening. However, the returned Appell FI values 
are not accurate for all input parameters, based on a comparison with Mathematica and 
the pure Python library mpmath. Even for cases where FI is correct, the computation is 
over an order of magnitude slower than numeric integration (for an error tolerance of 0.1 
ppm). I concluded that integration is a faster and easier solution than analytic models for 
the four-parameter nonlinear limb darkening law. 


4. Summary 

I introduced a new algorithm for computing transit light curves for any radially 
symmetric stellar limb darkening law. I also described the open-source Python package 
batman, a versatile code to generate model transit and eclipse light curves. Uniform, lin¬ 
ear, quadratic, logarithmic, exponential, and four-parameter nonlinear limb darkening laws 
are currently supported, batman uses C extension modules to compute light curves and is 
parallelized with OpenMP to optimize performance. Light curves can be calculated with 
accuracy better than 0.001 ppm, ensuring that the community is prepared to model the 
extraordinarily precise data we anticipate from upcoming facilities, batman is available at 
https : //github. com/lkreidberg/batman and is also hosted on the Python Package Index 
under the name batman-package. 


I thank Jacob Bean, Kevin Stevenson, Eric Agol, Ethan Kruse, Geert Jan Talens, 
Thomas Beatty, Brett Morris, and Karl Fogel for their support in developing batman. I 
also thank contributors to SciPy, Matplotlib, and the Python Programming Language for 
software and services. I am grateful for helpful suggestions from the referee, Jason Eastman, 
that improved the code and the manuscript. Support for this work was provided by a 
grant from the National Science Foundation through a Graduate Research Fellowship to the 
author. 
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